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ABSTRACT 

The leading hadronic contribution to the muon anomalous magnetic moment 
is given by a weighted integral over euclidean momentum of the hadronic vacuum 
polarization. This integral is dominated by momenta of order the muon mass. 
Since the finite volume in lattice QCD makes it difficult to compute the vacuum 
polarization at a large number of low momenta with high statistics (combined 
with the fact that one cannot compute it at zero momentum), a parametrization 
of the vacuum polarization is required to extrapolate the data. A much used 
functional form is based on vector meson dominance, which introduces model 
dependence into the lattice computation of the magnetic moment. Here we 
introduce a model-independent extrapolation method, and present a few first 
tests of this new method. 
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I. INTRODUCTION 



Recently, there has been an increased interest in the lattice computation of the hadronic 
contributions to the anomalous magnetic moment of the muon [1-6]. The aim is to 
provide a first-principles computation of the hadronic contributions to from lattice QCD 
assuming the Standard Model, with full control of the error. Since the experimental value of 

is known with great accuracy [7] , this would provide an interesting test of the Standard 
Model if the theoretical computation can be carried out with a comparable precision [8] . 

The dominant hadronic contribution comes from the hadronic vacuum polarization, and 
can be written as [1, 9]^ 

fiQ') 

where Ili{Q^), defined by 

n^u-AQ) = iQ\'^ - QMMQ') , (i-2) 

is the (flavor-diagonal) vacuum polarization for quark flavor i. In Eq. (l.ld), qi is the electric 
charge of quark i in units of the electron charge, a is the fine-structure constant, and 
is the muon mass. Disconnected contributions are expected to be small [10] and have been 
neglected in Eq. (l.ld).^ Of course, Eq. (1.1a) also applies to the electron and tau magnetic 
moments, if one replaces the muon mass by the electron or tau mass. 

In principle, the integral in Eq. (1.1a) could be approximated by a sum, if lattice values 
of n((5^) at sufficiently many low values of (of order m^, the region that dominates 
the integral) could be computed with high precision. At present, this is not yet possible. 
Instead, all lattice computations of a^^^ rely on fitting the lattice data to a functional form 
for n((5^), which is then used to compute the integral. The most successful approach to 
date has been to choose a functional form inspired by vector- meson dominance ( VMD) , with 
either only a contribution from the p pole (possibly dressed up with pion loop contributions 
[3]), or with contributions from the p and the p' [5], where the p mass is fixed to its value 
as independently determined on the same ensemble of gauge-field configurations. 

However, this introduces a model-dependent element in what is supposed to be a first- 
principles computation in lattice QCD. This results in a systematic error afflicting lattice 
computations of a^^° which is difficult to quantify.^ While VMD fits to Ii{Q'^) look very 



^ Our sign convention for n(<5^) is opposite to that of Ref. [3]. 

^ The methods developed in this article also apply if disconnected parts are included as well. For numerical 
c'stimates of disconnected contributions, see Ref. [4]. 

As already discussed in Ref. [3], chiral perturbation theory is of little help in this case. 



POO 

4a' J dQ'f{Q'){U{0)-U{Q')) , 

l + m2Q2^2(Q2) 

E ^^'n.(g^) , 

i=u,d,s 



(1.1a) 
(1.1b) 

(1.1c) 
(l.ld) 
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good, and lead to quite small statistical errors, is extremely sensitive to the behavior 
of the fitted Il{Q^) at very small Q^, and thus to any systematics affecting the small- 
behavior of Il{Q^). Therefore, it would be very nice if a functional form of Il{Q^) 
could be constructed that is solely based on known mathematical properties of the vacuum 
polarization, and which can be systematically improved if data with higher precision become 
available. 

It turns out that such a method exists. It is based on the well-known observation that 
the vacuum polarization can be expressed in terms of a positive spectral function through 
a (once-subtracted) dispersion relation. This makes it possible to express the vacuum po- 
larization n(Q^) in a form for which a convergent sequence of Pade approximants (PAs) is 
known to exist. Moreover, the convergence is uniform for any compact region in the com- 
plex plane excluding the cut along the negative real axis. This includes, in particular, any 
finite interval in euclidean between and oo. Since the contribution from say the region 
> 3 GeV^ to tt/j, is much smaller than currently attainable errors, this is sufficient to 
employ this observation for the computation of a^. Our goal is to present an exploration of 
this observation, using examples of available data for n(Q^). 

This article is organized as follows. In Sec. II we review elements of the necessary math- 
ematical theory, beginning with the observation that n(Q^) can be written in terms of a 
Stieltjes function, for which a converging sequence of PAs is known to exist. In Sec. Ill we 
explain how we will apply this construction in order to carry out fits to numerical data for 
Il{Q^). Section IV reports on two examples of such fits. We discuss future prospects of this 
approach in our concluding section. 

II. STIELTJES FUNCTIONS AND PADE APPROXIMANTS 

In this section, we review the necessary elements of the theory of Pade approximants (PAs) 
for functions that can be written as a Stieltjes integral with a finite radius of convergence. 
A good review is Ref. [11]; for multi-point PAs we refer to Refs. [12, 13]. 



A. Stieltjes functions 
Consider the function 



with z/(r) some real, bounded, non-decreasing function on the interval [0, 1/-R], taking in- 
finitely many values on that interval. The function ^{z) then is a Stieltjes function, and it 
is analytic everywhere in the complex plane except on the negative real axis for z < —R. 
The function ^{z = Q^) decreases monotonically as a function of for e {—R, oo). 

The vacuum polarization 11 (Q^) can be expressed in terms of $ through a once-subtracted 
dispersion relation 

n(Q2) = n(0) - Q'$(Q') , (2.2) 

Pit) 



$(g') = / dt 



4ml tit + Q') ' 
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where p(t) is the spectral function, which, of course, satisfies the constraint p(t) > for 
4m^ < t < oo. This can be seen by changing variables r = 1/t in the integral, taking 
R = 4m^, and choosing 

du{T) = p{l/T)dT , (2.3) 
p{t) = -lmn(t) . 

TT 

Let us consider an ordered sequence of positive values Qf of the variable Q^, with i G 
{1, . . . , P} and < Ql < Q2 < ■ ■ ■ < Qp, and assume that the function is known at 

these points. We may now construct a sequence of Stieltjes functions as follows. We begin 
by defining a function ^i{Q'^) by writing ^{Q^) as 

Then \E'i((5^) is also a Stieltjes function [12]. Moreover, "i/i^Q^) is positive on the interval 
[—R, 00), and, on that interval, has an upper bound 

< ^i(-R) < , Q' e [-i?, 00) . (2.5) 

This follows from the requirement that $((5^) not have a singularity on the real axis for 
> —R, which implies that 

lim {Ql - g')*i(g') < 1 . (2.6) 

Clearly, a sequence of Stieltjes functions ^^(Q^), i € {!,..., P}, can be constructed by 
iteration: 

= TT(^^ff|(Q5) . 'e{2.....P}. (2.7) 
which on the interval [— P, 00) satisfy 

<*.«^) <*.(-«) = ^ (1 - , .e{2.....P}, (2.8) 

where ^'i_i(-P) = limQ2^_^ ^i_i(g2). 

Equation 2.4 defines "ifiiQl) in terms of ^(Qi) and $(Qi)- Likewise, in general, ^i_i(Qi) 
can be expressed in terms of the values ^{Q^), j e {1, . . . , i} by using Eq. (2.7) recursively. 

Applying Eq. (2.7), the original function ^{Q^) can be written as a continued fraction 

1+(Q2-Q^)*p(Q2) 

As already observed above, ^'i(Qi), • • • , \l/p_i((5p) can be expressed in terms of the values 
of the function $(Q^) at the points Ql, . . . , Qp. 
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B. Multi-point Pade's 



A rational (or Padc) approximation to the function ^{Q^} can be constructed by setting 
^p(Q^) ill Eq. (2.9) equal to its lower bound {i.e., zero), or its upper bound, given by 
Eq. (2.8). A rational approximation R^{Q'^) is the ratio of two polynomials of degrees N 
and M 

We will refer to R^iQ^) as an [A^, M] PA.^ 

If we choose ^'i(Q^) = 0, the expression in Eq. (2.7) yields a [0,0] PA for '^i_i{Q^). 
Working back to the original function, this choice leads to a PA for the function $((5^)- If 
the number of points Qj, i E {1, . . . , P} is even, P = 2k, starting with ^p((5^) = yields a 
[k — l,k] PA. Indeed, for a [/c — 1, k] PA we need to solve for k coefficients a„ and k coefficients 
bn in Eq. (2.10), for a total of P = 2k, determined by the values ^{Qf), i G {!,••• ,P}- 
Likewise, for P = 2A; + 1 odd, the procedure yields a [k, k] PA. In short, from ^p{Q'^) — 
one obtains a [[(P — 1)/2J, [P/2J] PA, where [x\ denotes the integer part of x. 

These "standard" multi-point PAs were studied in Refs. [12, 13]. By construction, they 
arc exact at the vahics = Qf, i.e., the PA takes precisely the values ^{Qf) at these 
values of Q^. Moreover, these PAs converge to the function ^(Q^). More precisely, if we 
consider a sequence of standard multi-point PAs constructed from the values of ^{Q^) at 
a collection of points Qj < Ql < ■ ■ ■ < Qp < Ql < oo with hmp_^ooQp = Ql, the PAs 
converge uniformly on any closed and bounded domain in the complex Q^-plane excluding 
the cut -oo <Q'^ <-R [13], for P ^ oo. 

If we choose ^p{Q'^) equal to the upper bound of Eq. (2.8), the expression in Eq. (2.7) 
yields a [0,1] PA for \l'p_i((5^). Again working back to the original function, this choice 
also leads to a PA for the function $(Q^). Now if the number of values Qf, i G {1, . . . , P} is 
even, P = 2k, this yields a [k, k] PA. The counting argument is analogous to that above, but 
now this PA has, by construction, a pole at — —R, which provides the extra information 
needed to find the 2A;-|-1 coefficients in Eq. (2.10) for this case. Likewise, for P = 2A;-|-1 odd, 
the procedure yields a. [k,k + 1] PA. In short, in this case wc obtain a [[P/2J, [(P + l)/2j] 
PA. These "complementary" multi-point PAs are also exact at the values — Qf; they 
were introduced in Ref. [13]. 

If, given P values ^{Qf), i G {1, . . . , P}, the standard PA is written as 

defining the polynomials Ap{Q'^) and Bp{Q'^), the complementary PA can be written as [13] 

pLP/2j .^2^ ^ (P + Ql)Bp^,{-R)Apm + (g^ - Ql)Bp{-R)Ap_,m 
^L(P+i)/2j ^ {R + Ql)Bp^^{-R)Bp{Q^) + {Q^-Ql)Bp{-R)Bp^,{Q^) ' ^''''^ 

As already mentioned, both these PAs are exact at the points Qf. Moreover, the comple- 
mentary PAs have a pole at Q^ = —R, as can be seen from Eq. (2.12). Between the points 



* Redundancy between the coefficients a„ and is removed by choosing one of them equal to 1. Here we 
choose 6m = 1- 
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and Qf^i, as well as between — it! and Qf and between Qp and oo, the standard and 
complementary PAs provide an upper and lower bound to the original function $((5^).^ 
Which is the lower bound and which the upper bound alternates as one progresses through 
these P + 1 intervals from —R to oo [12, 13]. 



C. Parametrization 

All the poles of our standard PAs should have their poles on the negative real axis, at 
locations < —R. Indeed, one can prove [11, 13] that these PAs can be written in the 
form 

LP/2J 

n=l """^^ 

with Go = for P even, and 

a„>0, ne{l,...,LP/2j}, (2.14) 
bip/2i > bip/2\-i > ■ ■ ■ > bi > R . 

Once these parameters have been obtained for the [|_(P — 2)/2j, |_(P — 1)/2J] and [|_(P — 
1)/2J,LP/2J] PAs, the complementary [L-P/2J, [(-P + l)/2j] PA can be obtained from 
Eq. (2.12). 



III. FIT STRATEGY 

In the situation of an actual fit to values of Il{Q^) obtained from a numerical computation, 
these values are only known within some statistical errors. That implies that we do not know 
any points of the function exactly. However, the fact remains that n(Q^) can be expressed in 
terms of a Stieltjes function, cf. Eq. (2.2), and this implies that a scries of PA representations 
as described in Sec. II B exists that converges to Il{Q^), when we parametrize the vacuum 
polarization as 

U{Q^) = n(0) - (^ao + f2 . (3.1) 

A number of PAs can be estimated by fitting this form to the data as a function of increasing 
N. For Oo = the parameters to be fitted are 11(0) and the a„ and 6„ for n e {1, . . . , N}, 
and we obtain an [iV — 1, A^] PA. When also is fitted we obtain an [A^, A^] PA, and together 
this sequence of PAs estimates the standard PAs introduced in Sec. II B. We will enforce 
the restrictions (2.14) on the parameters a„ and 6„ in our fits (unless stated otherwise). 

Of course, since these PAs will not have been constructed as exact multi-point PAs, the 
bounds described at the end of Sec. II B will not be exact either. However, we can still check 
the stability of our results as a function of A^. In particular, we can check the stability of the 
value we obtain as a function of A^ for the quantity a^^*^ defined in Eq. (1.1a). As we will 
see, this quantity is much more stable as a function of A^ than the individual fit parameters 



^ These bounds are optimal [12, 13]. 
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X /dot 


inlO HLO,Q^<l 
ID - 


n(0) 




h 




VMD 


5.86/3 


363(7) 


0.0962(6) 


0.0471(9) 


0.9256(nxed) 




[0,1] 


11.4/8 


338(6) 


0.0960(5) 


0.0600(7) 


1.287(27) 




[1,1] 


7.49/7 


350(8) 


0.0963(6) 


0.049(4) 


1.09(9) 


0.0028(12) 


[1,2] 


7.49/6 


350(8) 


0.0963(6) 


0.049(4) 


1.09(9) 












2(17) 


2(8) X 10^ 




[2,2] 


7.49/5 


350(7) 


0.0963(6) 


0.049(4) 


1.09(9) 


0.0012(10) 










2.4(1.4) 


1.4(0.8) X 10^ 





TABLE 1: VMD and PA fits to the a = 0.09 fm, amught = 0.0124 data for U{Q'^) of Ref. [3] with 
Q"^ < 0.6 GeV^, except for the VMD fit, for which the fit interval is < 0.35 GeV^. Correlated 
fits; errors. 

in Eq. (3.1). In practice, we have explored fits up to A?" = 3; here we will present fits up to 
N = 2. 

As in Ref. [3], we will mostly explore fits in which we take the values of for which we 
fit the PAs from the data for Yl{Q'^) in an interval between — and — 1 GeV^. For 
each of our fits, we will compute the quantity 

pl GeV^ 

^HLo,Q^<i ^ / f^Q2^ (n(o) - n(g^)) , (3.2) 

^0 

with f{Q'^) defined in Eq. (1.1b). This of course misses the part of the integral between 
1 GeV^ and oo, but this part is of order a percent of the low-Q^ contribution. Since our 
goal here is to test the Pade approach to fitting Il{Q^) for < 1 GeV^, we have restricted 
ourselves to the expression in Eq. (3.2) for comparisons between different fits. 

We have not explored the complementary PAs defined in Sec. 11 B yet, but we anticipate 
that they may become useful in the future, when more precise data become available. 

IV. TESTS 

In this section, we explore fits to two different data sets. One set is the data for Il{Q^) 
with hght quark mass equal to 0.0124 in lattice units (on a 24^ x 96 lattice with lattice 
spacing a ~ 0.09 fm, using "fine" configurations from the MILC collaboration [14]) that 
was also studied in Ref. [3] (see Table 1 of that paper). The other set is data obtained using 
the MILC "super-fine" gauge configurations with lattice spacing a ~ 0.06 fm on a 64^ x 144 
lattice with light quark mass equal to 0.0018 and a strange quark mass equal to 0.018, in 
lattice units. For both data sets the lattice strange quark mass is approximately equal to 
the physical strange quark mass. We will always assume that the theory of Sec. 11 applies to 
these data, i.e., that lattice artifacts are small enough to be ignored. In practice, rotational 
invariance is broken on the lattice. Since Tl{Q^) is extracted from Eq. (1.2) using a lattice 
definition of the momentum components [3], breaking of rotational invariance causes 
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X /dot 


inlO HLO,Q^<l 
ID - 


n(o) 




L 

k 


ao 


VMD 


4.37/18 


413(8) 


0.0980(7) 


0.0536(10) 


0.9256(nxed) 




[0,1] 


3.58/17 


373(37) 


0.0971(12) 


0.0569(25) 


1.10(16) 




[1,1] 


3.36/16 


424(116) 


0.0979(22) 


0.033(14) 


0.6(4) 


0.007(6) 


[1,2] 


3.35/15 


443(293) 


0.098(4) 


0.02(10) 


0.4(1.7) 












0.058(12) 


2(11) 




[2,2] 


3.35/14 


445(432) 


0.098(4) 


0.02(29) 


0.4(4.2) 


0.0(4) 










0.1(3.8) 


4(141) 





TABLE 2: VMD and PA fits to the a = 0.09 fm, amught = 0.0124 data for n(Q2) of Ref. [3] with 
< 1 GeV^. Uncorrelated fits 

. QjYors computed by a linear fluctuation analysis. 

n((5^) to show small deviations from the monotonic decrease that follows from Eq. (2.2).^ 

For each data set we carry out both correlated and uncorrelated fits, and compare those 
with each other. It turns out that this raises interesting questions about the behavior of the 
data and the fits at very low Q'^. 

A. a — 0.09 fm data at mught/m strange — 0.4 

For our first example we consider the amnght = 0.0124 data that were also considered 
in Ref. [3]; this value of the light quark mass corresponds to about 2/5 times the physical 
strange quark mass. In Tables 1 and 2 we show the result of a sequence of PA fits, with 
Tabic 1 showing correlated fits, and Tabic 2 showing uncorrelated fits. For the correlated fits 
we fitted data on the interval < < 0.6 GeV^, because this interval yields the smallest 
values for the V^^ degree of freedom. For the uncorrelated case, we fitted the data for 
n((52) on the interval < < 1 GeV^, as was done in Ref. [3]. 

Table 1 shows that the value for a^^'^''^ becomes very stable for PA fits starting at 
[1,1]. For the [1,1] PA x^/dof = 1.07, indicating a good fit. For higher PAs, the value 
of does not change, being very insensitive to the location of the second pole. It follows 
that the values of the parameters characterizing the second pole are not well determined, as 
can be seen in the table. However, the value of a^^'^''^ is completely insensitive to the 

second pole. The explanation for this is that the integral for a^^^''^^-^ is dominated by the 
Q'^ region around m^, and thus very insensitive to the precise location of PA poles at large 

negative values of Q^. 

The fit marked "VMD" is obtained by holding the parameter bi fixed at the square of the 
p mass (which is equal to 962 MeV for this data set) in what would otherwise be a [0, 1] PA 
fit. It is thus not one of the sequence of PAs introduced in Sec. II. According to the theory 
there is no reason one should expect the parameter bi to be equal to the square of the p 
mass, as borne out by the values for bi found in the PA fits of Table iJ For the correlated 



^ Also, all data have been obtained on cnsemblcis of gauge configurations generated with improved actions. 
Therefore, the spectral function may not be positive for values of the momenta near the lattice cutoff. 
See Sec. IV C for further discussion. 
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0.100 




Q' (GeV' ) 

FIG. 1: [1, 1] fits of Tables 1 (correlated, solid curve) and 2 (uncorrelated, dashed curve) compared 
with data. Solid points have been included in the correlated fit while both solid and open points have 
been included in the uncorrelated fit. 

VMD fit a fitting interval < < 0.35 GeV^ leads to the lowest per degree of freedom. 
With x^/dof ~ 2, the VMD fit is not very good. It is already much better for the [0, 1] PA, 
in which the constraint on hi is relaxed, and it decreases further, to an acceptable value, for 
the [1, 1] PA. 

Table 2 shows similar fits, but here all fits are uncorrelated. All errors have been esti- 
mated using a linear fluctuation analysis starting from the uncorrelated y^, starting from 
the full data covariance matrix [15]. These errors agree with errors computed under a single- 
elimination jackknife. In these PA fits we have relaxed the constraint hi > 4m^ = 0.906 GeV^ 
(on this data set), but one notes that the values of hi are consistent with this bound within 
errors. Both correlated and uncorrelated [1, 1] PA fits are shown in Fig. 1. 

The uncorrelated VMD fit reproduces "fit A" of Ref. [3], including the error. ^ One notes 
that the uncorrelated PA fits lead to results consistent with those of Table 1, but with much 
larger errors. The uncorrelated VMD fit is not consistent with what we would expect to be 
the best fit, 

from the [1, 1] PA of Table 1. 

We may also compare the values in the tables with values obtained from a fit with a 
fourth order polynomial in Q^, which are 

qHlo,q2<i _ 410(91) X 10-^° , (uncorrelated) , (4.2) 
HLO,Q2<i ^ 34g^g^ ^ ^Q-io ^ (correlated) . 



The parameters 11(0) and ai are not the same as the parameters A and fv of Ref. [3]. 
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X /dot 


inlO HLO,Q^<l 
ID - 


n(0) 




k 


ao 


VMD 


38.6/18 


646(8) 


0.1222(6) 


0.0595(8) 


0.64 (fixed) 




fn -I 1 

[0,1] 


14.3/17 


550(20) 


0.1203(7) 


0.0646(16) 


0.83(5) 




r-i -1 1 

[1,1] 


13.9/16 


572(41) 


0.1206(8) 


0.052(16) 


0.68(20) 


0.005(7) 


[1,2] 


13.9/15 


572(37) 


0.1206(8) 


0.052(14) 


0.68(19) 












1(6) 


0.3(1.0) X 10^ 




[2,2] 


13.9/14 


572(38) 


0.1206(8) 


0.052(14) 


0.68(18) 


0.003(27) 










1(31) 


0.4(6.0) X 10^ 





TABLE 3: PA fits to the a = 0.06 fm, amught = 0.0018 data for U{Q'^) with < 0.53 GeV^. 
Correlated fits; errors. 





xVdof 


^010^HLO,Q^<1 


n(0) 


ai 


bi 


ao 


VMD 


37.2/51 


685.2(7.8) 


0.1236(6) 


0.0631(7) 


0.64 (fixed) 




[0,1] 


13.9/50 


555(22) 


0.1208(8) 


0.0666(7) 


0.85(4) 




[1,1] 


12.0/49 


645(66) 


0.1221(13) 


0.047(5) 


0.54(11) 


0.0071(21) 


[1,2] 


11.4/48 


788(482) 


0.123(4) 


0.015(20) 


0.2(4) 












0.063(14) 


1.4(9) 




[2,2] 


11.3/47 


837(627) 


0.124(5) 


0.018(5) 


0.2(5) 


0.022(9) 










0.22(6) 


3.9(6) 





TABLE 4: PA fits to the a = 0.06 fm, amught = 0.0018 data for U{Q^) with < 1 GeV^. 
Uncorrelated fits; errors from linear fluctuation analysis. For the [1,2] and [2,2] fits, hi is at the 
limit Arn^ = 0.1936 GeV^ (for this ensemble), which was enforced in those fits. 

The first line is in agreement with Ref. [3], and was fitted with < < 1 GeV^, as in 
Table 2, and the second is from a correlated fit on the interval < < 0.6 GeV^, as in 
Table 1. The latter fit has a x^/dof of 7.48/6, less good than the [1, 1] fit in Table 1. Both 
are in good agreement with Eq. (4.1), given the size of the errors. 

B. a — 0.06 fm data at mught/m strange — 0.1 

For our second example, we consider the vacuum polarization computed on MILC con- 
figurations at a = 0.06 fm and amught = 0.0018, which is about 1/10 times the physical 
strange quark mass. Correlated fits are shown in Table 3, where we fitted the data for 
< Q'^ < 0.53 GeV^ (which corresponds to the 20 data points with the lowest values of 
Q^). The values per degree of freedom of the fits in Table 3 arc slightly smaller than 
one, except for the VMD fit, for which x^/dof is about two.^ We find that the value of 
X^/dof increases if we fit over a larger range of values, and we will therefore take the 



We thank Doug Toussaint for providing us with an unpublished rough estimate of the p mass for this data 
set. 
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0.13 




0.08' ' ' ' ' 1 

0.0 0.2 0.4 0.6 0.8 1.0 

(GeV' ) 



FIG. 2: [1,1] PA fits of Tables 3 (correlated, solid curve) and 4 (uncorrelated, dashed curve) 
compared with data. Solid points have been included in the correlated fit while both solid and open 
points have been included in the uncorrelated fit. 



results ol Table 3 as our optimal results (for more on this point, see the discussion around 
Table 5 below). Uncorrelated fits are shown in Table 4, where, in line with Sec. IV A, fits 
were carried out on the interval < < 1 GeV^. 

It is again not surprising that the correlated fits become less good if one fits over a larger 
range in Q^. As before, it is clear from the tables that, given the quality of the data, it is 
very hard to fit a second pole. The value of a}^^^'^ is again completely insensitive to the 
location of the second pole.^° 

We show the [1, 1] fits of Tables 3 and 4 in Fig. 2. As in Fig. 1 one notes the sensitivity 
of the fit near = 0; this explains the different values for aJJ^^'*^ shown in the tables. 

From the [1, 1] PA fit of Table 3 we take what we would expect to be our best result for 
this data set: 

^hlo,q2<i ^ 572(41) X 10"^° . (4.3) 

In Fig. 3 we show correlated and uncorrelated [1, 1] PA fits, now taking the range < < 
0.53 GeV^ as our fitting range also for the uncorrelated fit. We note that the uncorrelated 
fit appears to do better than the uncorrelated [1, 1] PA fit shown in Fig. 2 at the lowest 
value, but much less well than the correlated fit for > 0.53 GeV^. Accordingly, 
uncorrelated fits are quite sensitive to the fitting range. For instance, the central value of 
^HLO,Q <i fj-QYa the uncorrelated fit shown in Fig. 3 is 42% larger than from a similar fit on 
the range < < 1 GeV^ (shown in Table 4). A correlated fit on the latter range gives a 
central value which is only 3% larger than the value in Eq. (4.3), i.e., it is within the error 



We even considered [2, 3] and [3, 3] fits, with the conclusion being the same. 
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FIG. 3: [1,1] correlated (solid curve) and uncorrelated (dashed curve) fits as in Table 3, fitted on 
interval < < 0.53 GeV^. Solid points have been included in the fits, open points have not been 
included. 



given in that equation. 

In Ref. [3] also polynomial fits with third- and fourth-order polynomials were considered, 
and it is thus interesting to compare PA fits with polynomial fits. For the data of this 
subsection, the radius of convergence, 4m^ = 0.194 GeV^.^^ Therefore, fitted polynomials 
cannot be interpreted as estimates of the Taylor expansion of Il{Q^) around = 0, as long 
as we use a fitting interval with upper bound larger than 4m^. 

We show third- and fourth-order polynomial fits, as well as [1,1] and [1,2] PA fits in 
Table 5, as a function of the number of data points in the fit (20 points corresponds to 
the fitting interval < < 0.53 GeV^ used in Table 3). All fits shown are correlated. 
Both "Poly 3" and "PA [1,1]" are four-parameter fits, while "Poly 4" and "PA [1,2]" are 
five-parameter fits. The x^/dof for all fits is good, except for fits with 26 data points, for 
which it shows a steep increase. 

We observe that Poly 3, PA [1,1] and PA [1,2] fits all lead to values for aHLO,Q2<i ^-^ich 
are stable within the error given in Eq. (4.3). For the Poly 4 fit, however, this spread is 
much larger. Adding a fit parameter by going from Poly 3 to Poly 4 fits leads to significant 
changes in the central value for a^'"^''^ whereas going from PA [1,1] to PA [1,2] fits the 
central values do not change much. In other words, if we would do a correlated Poly 4 fit 
to 20 data points, for which we would find a^^'^''^ = 535(45) x 10"^*^, the error would be 



Despite the fact that for a correlated fit on the range < Q'^ < 1 GeV^ the value of is about 2.5 per 
degree of freedom. 

The radius of convergence for the case of Sec. IV A is much larger, which is why we chose to make this 
comparison in this subsection. 
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Poly 3 


Poly 4 


PA [1,1] 


PA [1,2] 


# points 


x7dof 


(1) 

ay 


xVdof 


(1) 

ay 


9/1 n 

xVdof 


(1) 

ay 


9/1 f 

xVdof 


(1) 

ay 


16 


9.6/12 


543 


9.5/11 


483 


9.7/12 


564 


9.7/11 


565 


18 


11.4/14 


526 


10.5/13 


596 


11.2/14 


541 


11.5/13 


561 


20 


13.1/16 


536 


13.1/15 


535 


13.9/16 


572 


13.9/15 


572 


22 


16.5/18 


541 


15.9/17 


513 


18.5/18 


566 


18.5/17 


566 


24 


16.6/20 


537 


16.4/19 


521 


19.4/20 


583 


19.4/19 


583 


26 


30.7/22 


505 


23.6/21 


580 


26.8/22 


557 


26.7/21 


560 



TABLE 5: Correlated PA and polynomial fits to the data of Table 3, as a function of the fitting 
interval. The first column shows the number of data points in the fit, with 20 points corresponding 
to the fitting interval < < 0.53 GeV^ of Table 3. The column "Poly n" shows results from a 
fit to a polynomial of degree n; a^^^ stands for 10^''a|^^^''^ . 
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FIG. 4: Comparison of correlated [1,2] PA (solid curve) and 4th- order polynomial (dashed curve) 
fits, both fitted on the interval < < 0.53 GeV^. Solid points have been included in the fits, 
open points have not been included. 



underestimated because of the spread of values for the Poly 4 fit shown in Table 5, while 
the error shown in Eq. (4.3) encompasses the full range of a^^'~'''^ PA results shown in 
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Table 5. The Poly 4 and PA [1,2] fits with 20 data points are shown in Fig. 4. 



C. Discussion of fits 

In this subsection we will discuss the fit results presented in Tables 1 through 4 in more 
detail. We begin with the a = 0.09 fm results of Tables 1 and 2. 

It is important to emphasize again that the VMD fits are not part of the sequence of 
PAs introduced in Sec. 11, because in the VMD fits the pole at = —hi is held fixed at 
the estimated (squared) p mass on this ensemble. The actual QCD spectral function has a 
cut on the negative axis starting at Q'^ — — 4m^; any poles reside on the second Riemann 
sheet, away from the negative axis.^^ While positivity of the spectral function imphes that 
n(Q^) can be expressed in terms of a Stieltjes function, with a convergent sequence of PAs 
given by Eq. (2.13), there is no reason that any of the poles of these Pade's should be equal 
to (the real part of) any pole representing a resonance in QCD. In particular, in the PA 
fits, the parameter hi should not be taken equal to the square of the p mass, but instead it 
should be left as a free parameter. We included the VMD fits in Tables 1 and 2 in order to 
compare them with the PA fits. 

First, we note that the correlated VMD fit in Table 1 is a rather poor fit, with a high 
X^/dof, and there is no agreement between the correlated and uncorrelated VMD fits. The 
quality of the correlated fits improves when we add more parameters, first by varying hi and 
then by adding in the parameter by which time x^/dof 1. 

The next observation is that the correlated fits do not get better by adding a second pole 
to the PA. The minimum value of stays the same, and consequently, the parameters of 
the second pole are very poorly determined. We have checked that this does not depend on 
the fitting range employed. A possible explanation is that the lattice data for the vacuum 
polarization do not quite follow the behavior predicted by Eq. (2.2), because of the breaking 
of rotational invariance on the lattice (c/. discussion at the beginning of this section). This 
can be seen from the fact that the data points in Fig. 1 show small deviations from a 
smooth, monotonically-decreasing behavior. However, we note that the value of a^^'^ '^ 
is completely insensitive to the parameters of the second (and higher) poles. Our best 
correlated fit value for a)^^^''^ is given in Eq. (4.1). 

Moving to the uncorrelated fits of Table 2, we observe that the values for a^^^''^^-^ for 
all fits in that table are consistent with each other, because of the rapid increase of errors 
with the order of the PA. Furthermore, all uncorrelated and correlated PA fits are consistent 
with each other as well, but clearly the correlated fits have much smaller errors. 

It is instructive to compare the best correlated fit, the [1, 1] PA fit, with the uncorrelated 
VMD fit, because both have very small errors, and provide a good fit, as can be seen in 
Fig. 5. Both fits are good fits, but they lead to values for a^^'^'^ which are not consistent 
with each other. The statistical error on the uncorrelated VMD fit is very small, but this fit 
has an unknown systematic error because of its model dependence. This may explain the 
discrepancy with the correlated [1, 1] PA fit. The latter may be expected to have a much 
smaller systematic error, since it is known to be a member of a converging sequence of PAs, 
and clearly already provides a good fit also in the larger region that was not included in 



Only in the limit of an infinite number of colors do the poles move toward the negative real axis, and the 
vacuum polarization becomes a meromorphic function. 
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0.100 




Q' (GeV' ) 

FIG. 5: [1,1] fit of Table 1 (correlated, solid curve) and VMD fit of Table 2 (uncorrelated, dashed 
curve) compared with data. Solid points have been included in the correlated fit while both solid and 
open points have been included in the uncorrelated fit. 

the fit (the horizontal axis of Fig. 5 covers about twice the fitted region < < 0.53 GeV^). 

However, for the computation of o}^^'^'^ the region ~ m^^ = 0.011 GeV^ dominates, 
and it is clear that the data do not distinguish between these two fits in that region. While 
the rapid increase of the goodness of fit seen in Table 1 from the VMD fit to the [1, 1] fit can 
be taken as an indication that correlated fits unbiased by model dependence are promising, 
we conclude that it is not possible to exclude either value of the basis of these 

data. 

Similar remarks apply to the fits shown in Tables 3 and 4. Correlated PA fits all have 
X^/dof ~ 1, unlike the correlated VMD fit for which x^/dof ~ 2. In both cases, the 
correlated and uncorrelated VMD fits do not agree within errors (which, we recall, are purely 
statistical). The value we obtain for 

^HLO,Q^<l ig 

larger than that obtained in Sec. IV A; we 
believe that this is mostly due to a smaller pion mass, with itIt^ ~ 220 MeV for this data 
set, while m^r ~ 480 MeV for the a = 0.09 fm data set. 

V. CONCLUSION 

In this article, we presented a new, model-independent method for fitting the hadronic 
vacuum polarization Il{Q^) as a function of euclidean momentum to data obtained from 
a lattice QCD computation. The method is based on the theory of Fade approximants (PAs) 
to a Stieltjes function, and yields, in principle, a converging sequence of FAs to the vacuum 
polarization. 

These FAs can be used to obtain lattice estimates for the leading hadronic contribution to 
the anomalous magnetic moment of the muon from Eq. (1.1a). By comparing successive FAs 
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in the sequence, one should be able to check the convergence in practice. This would allow 
for a fully model- independent determination of the leading hadronic contribution a^^*^, and 
thus help eliminate an unknown systematic error present in all lattice computations of a^^'^ 
to date. 

In comparison with the VMD fits which have been employed in the literature, these 

PAs contain more parameters (the [0, 1] PA already contains three parameters, whereas the 
simplest VMD ansatz contains only two). One thus typically expects larger statistical errors 
given certain lattice data. However, the PA approach avoids model-dependent assumptions, 
and hence removes the unknown systematic error associated with the VMD approach. 

We have explored this new framework on two state-of-the-art ensembles of gauge configu- 
rations, at different lattice spacings and pion masses. From these explorations, we conclude 
that this new method looks promising, but that better data at very low values of will 
be needed in order to control the extrapolation necessary for a reliable computation of aJJ'"*^ 
from the integral in Eq. (1.1a). 

Our explorations show that given current lattice data for n(Q^), there is a significant 
difference between our best PA fits (which are four-parameter [1, 1] PAs), and VMD fits.^^ 
For instance, the difference between the values of aJJ^*^'*^ obtained from the correlated 
[1, 1] PA fit of Table 1 and the uncorrected VMD fit of Table 2 and Ref. [3] is about 15-20%, 
much larger than the statistical fit error on each of these values. While it is tempting to 
view the value from the correlated [1, 1] PA fit as the more reliable one, it is clear from 
Fig. 5 that more data points with higher precision at low are needed in order to reduce 
this uncertainty. 

Our explorations also showed that with these data it is very difficult to fit the parameters 
characterizing the second and higher poles of the PAs. In order to test the convergence of 
the sequence of PAs fitted to II (Q^), it would be desirable to investigate this issue, which 
is possibly related to breaking of rotational invariance at non-zero lattice spacing, in more 
detail in the future. This issue appears to have no direct effect on the value of a}^^'^'^ 
which we found to be very insensitive to the location and residues of the second and higher 
poles. 

In conclusion, the new method presented here looks promising, but data for Il{Q^) with 
more values at ~ and with higher statistics will be necessary in order to attain the 
high precision determination of a^'^ needed for a meaningful comparison with experiment. 
Work in this direction is in progress. 
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